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Abstract 

The dynamical effects of the tensor modes of the geometry are investigated in the context 
of curvature bounces. Since the bouncing behaviour implies sharp deviations from a radiation- 
dominated evolution, significant back-reaction effects of relic gravitons may be expected at short 
wavelengths. After developing a general iterative framework for the calculation of dynamical 
back-reaction effects, explicit analytical and numerical examples are investigated for different 
parametrizations of the energy-momentum pseudo-tensor(s) of the produced gravitons. The 
reported results suggest that dynamical back-reaction effects are a necessary ingredient for a 
consistent description of bouncing models at late times. 
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1 Introduction 



Relic gravitons are copiously produced in the early Universe due to the pumping action of the back- 
ground geometry [1,2]. If a quasi-de Sitter phase of expansion is followed by a radiation-dominated 
phase, the logarithmic energy spectrum (in units of the critical energy density), customarily de- 
noted by Q gw (u), is quasi-flat [3] (see also [4, 5, 6]) for present (physical) frequencies v ranging 
between 1CP 16 Hz and, approximately, 100 MHz. The transition from the radiation-dominated to 
the matter epoch produces an infrared branch where O gw ~ v~ ~ 2 between 10" 18 Hz and 10~ 16 Hz 
[7, 8, 9]. 

In pre-big bang models [10, 11], the spectrum of relic gravitons is far from scale-invariant. The 
spectral slope of Q gw (v) (for frequencies larger than 10~ 16 Hz) is violet, i.e. Q gw ~ v" f with 7 > 1 
[5, 12] (see also [13] and, in a complementary perspective, [14]). Minimal pre-big bang models are 
characterized by a slope 7 = 3 that become progressively less steep as the frequency increases in 
the kHz and MHz region. 

Depending on the features of the cosmological model, the energy density of relic gravitons 
may affect the dynamics and change the time evolution of the scale factor as well as of the other 
homogeneous quantities. If the equation of state of the sources of the background geometry is 
stiffer than radiation (i.e. p = wp with w > 1/3) the relativistic gravitons of short wavelengths can 
change the dynamical evolution, as it was noticed in the context of the stiff model of Zeldovich [2]. 

When a stiff phase follows an inflationary stage of expansion [15] (see also [16, 17] 2 ) the 
back-reaction effects of the produced gravitons set a limit on the possible duration of the stiff 
(post-inflationary) phase. This aspect is relevant in different situations and, for instance, in the 
context of quintessential inflationary models [20] where Q gw oc v for frequencies larger than the mHz 
[22, 23, 24] (see also [25, 26, 27, 28] for possible detection strategies). Blue spectral slopes have been 
also derived in the context of quintessential models on the brane [29, 30] (see also [31]). Constraints 
arising from the production of relic gravitons in quintessential models have been discussed in [32]. 

Defining as r the conformal time coordinate and k as the comoving wave-number, two physically 

different regimes appear naturally in the problem. If kr <C 1 (short wavelength limit) the Fourier 

modes of the tensor fluctuations of the geometry are said to be super-adiabatically amplified. In 

the opposite regime, i.e. kr S> 1 (long wavelength limit) the tensor modes are oscillating. In the 

2 It is relevant to remark that the authors of Refs. [16] and [17] choose to compute the energy and pressure 
density of the produced particles by means of a perturbative expansion whose small parameter is the deviation from 
conformal coupling (|£ — 1/6| in the notations of [17]). This approach (also discussed in [18, 19]) may be applied, with 
some caveats, in the minimally coupled case (corresponding to £ — » 0) where, it could be argued that the expansion 
parameter does not exceed 1. In the present investigation this expansion will not be used since the coupling of relic 
gravitons to the background geometry is not close to conformal. Furthermore the expansion in |£ — 1/6| may partially 
break down during the bouncing regime. When needed, numerical techniques will be employed. 
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short wavelength limit the effective equation of state obeyed by the relic gravitons is radiative, i.e. 
(Pgw) = (/0gw)/3 where (p gw ) and (pgw) correspond to the averaged pressure and energy densities. 
In this case, (p gw ) ~ a -4 , where a(r) is the scale factor. In the long wavelength limit the effective 
equation of state may depend on the salient features of the background evolution. For instance, 
during the transition from a quasi-de Sitter stage of expansion to a radiation dominated stage, the 
modes of long wavelengths are compatible with an effective equation of state (p gw ) = — (Pgw)/3 
implying that (p gw ) ~ a~ 2 . In this limit the system behaves, effectively, as a generalized fluid 
dominated by the spatial gradients of the tensor modes of the geometry. 

While the conclusions mentioned in the previous paragraph must rely on specific definitions 
of the energy and pressure densities of the relic gravitons, it is well known that, in general, it is 
impossible to assign a localized energy density to the gravitational field [33]. This caveat does 
not exclude the possibility of adopting consistent frameworks for the analysis of a gravitational 
energy-momentum pseudo-tensor. 

To be more specific, the tensor modes of the geometry can be characterized by a rank-two tensor 
hij defined in the three spatial dimensions (that will taken to be flat), i.e. 

ds 2 = a 2 (r)[dT 2 -(Sij + hifidjdaP], h\ = dih) = 0, (1.1) 

where Sij is Kroeneker symbol and hij, being traceless and divergence- less, carries two independent 
degrees of freedom that correspond to the two polarizations of a gravitational wave in a conformally 
flat Friedmann- Robertson- Walker (FRW) background 3 . The simplest way of estimating the impact 
of the created gravitons on the background dynamics is by computing the lowest-order nonlinear 
corrections to the Einstein tensor 

q; = k- I^r, (i.2) 

where R 1 ^ and R are, respectively, the Ricci tensor and the Ricci scalar. The nonlinear corrections 
to the Einstein tensor, will consist, to lowest order, of quadratic combinations of hij that can be 
formally expressed as 4 

4t; = -j t (2) ^, (i.3) 

where the superscript at the right hand side denotes the second-order fluctuation of the correspond- 
ing quantity while the subscript refers to the tensor nature of the fluctuations. This procedure is 
essentially the one described in [34, 35] and has been re-explored, in a related context, by the 

3 Since hij is transverse and traceless, the direction of propagation can be chosen to lie along the third axis and, 
in this case the two physical polarizations of the graviton will be h\ = —h\ = /i® and h\ = h\ = h®. This will be 
the nomenclature followed in the present paper. 

4 In the present notations the Planck length will be defined as £p = \/8ivG in units h — c = 1. Natural gravitational 
units 16itG = 2l\ — 1 will also be adopted when needed. 
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authors of Refs. [36, 37] mainly in connection with conventional inflationary models where the 
Universe is always expanding 5 

In [39] (see also [40]) a complementary perspective was invoked and it was observed that, by 
perturbing the Einstein-Hilbert action to second order in the amplitude of the tensor modes of the 
geometry in a Friedmann-Robertson- Walker background, each single polarization of the graviton 
obeys, up to total derivatives, the action of a minimally coupled scalar degree of freedom 

5l 2) S =\J d A x a 2 (r) rf^d a hd p h, (1.4) 
where the dimensionless amplitude h has been defined as 

h= K = K_ (L5) 

Consequently, following [39] , a natural ansatz for describing the back-reaction of the created gravi- 
tons on the background space-time would be to use, for each single polarization, the energy- 
momentum tensor of a minimally coupled scalar field. This approach has been exploited, for 
instance, in [6] and in [17] for estimating the energy and pressure densities of the relic gravitons 
in a multi-stage Universe. A possible weakness of this second approach is that Eq. (1.4) holds, 
strictly speaking, in a Friedmann-Roberston- Walker background. In spite of this, it is clearly useful 
to bear in mind this possibility. 

The purpose of the present paper is to scrutinize the dynamical back-reaction effects of the relic 
gravitons in the extreme situation where the background geometry undergoes a contraction that 
smoothly evolves, through a bounce, into an expanding phase. This sequence of events takes place, 
for instance, in pre big-bang models. The bouncing beahaviour is related to strong deviations from 
a radiative Universe. Consequently, short wave fluctuations are expected to modify or even restrict 
the dynamical range of bouncing models since, in this regime, the effective equation of state of relic 
gravitons is the one of radiation. For long wavelength fluctuations, the qualitative expectation is 
more difficult to formulate since, as already mentioned, the effective equation of state is sensitive 
to the specific features of the background evolution before and after the bounce. 

Consider then, as an example, one of the simplest realizations of bouncing dynamics, i.e. the one 

provided by minimal pre-big bang models where, in the Einstein frame, the Universe first contracts 

as a(r) ~ sf—rfrx and then expands, in the post-big bang phase, as a(r) ~ \J~tJt\- An accurate 

description of the bouncing regime in the region — n < r < t\ may be achieved, for instance, by 

means of a non-local potentials [41, 42, 43] that depend on the so-called shifted dilaton which is 

5 Recently in [38] the field theoretical formulation of General Relativity was further developed with the purpose of 
deriving an energy-momentum tensor for the gravitational field. Our approach, in the present context, is more modest 
since we simply want to compare the dynamical consequences of different possible definitions of the energy-momentum 
pseudo-tensor of the relic gravitons. 
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usually introduced in the context of the scale-factor-duality symmetry [10] (or more generally, in 
the treatment of 0(d,d) transformations [44, 45] acting on the background fields of the low-energy 
string effective action). The features of the regular bouncing solutions discussed in [46, 47] are such 
that the role of the dilaton potential is only crucial in the bouncing region while, away from the 
bounce (i.e. |r| > t\) the the geometry is driven by the kinetic energy of the dilaton field. 

As soon as the tensor modes reach into the super-adiabatic regime, i.e. kr <C 1 the effective 
pressure density of the relic gravitons becomes (p gw ) c± (p gw )- Consequently, at least up to the 
bounce, the energy density of the relic gravitons is not likely dominate on the energy density of the 
dilatonic sources since the two components evolve exactly in the same way with the scale factor. 
While it is plausible that only around the bounce the short wavelength modes become dynamically 
relevant, it is also clear that the whole qualitative picture should be corroborated by a detailed 
numerical and analytical treatment. 

It will be instructive to conduct the calculation of the dynamical back-reaction within different 
ansatz for the energy-momentum pseudo-tensor. Then the results will be compared. A byproduct 
of the present analysis will indeed be that different forms of the energy-momentum pseudo-tensor 
lead, quantitatively, to the same back-reaction effects. 

The present paper is then organized as follows. In Section 2 a class of string inspired bouncing 
cosmologies will be introduced. Section 3 is devoted to the analysis of the production of relic 
gravitons and to the different definitions of their energy-momentum pseudo-tensor in the framework 
of bouncing solutions. Section 4 deals with the analytical estimates of the effective barotropic 
indices. Section 5 describes the implementation of a self-consistent (iterative) scheme that allows 
to compute the back-reaction effects of relic gravitons. Section 6 summarizes the main lessons to 
be drawn from the present analysis and contains also the concluding remarks. Finally, various 
technical results that are relevant for the present investigation are reported in the appendix. 

2 A class of regular curvature bounces 

Regular bouncing solutions can be obtained by means of dilaton potentials that respect the scale 
factor duality symmetry [10]. These potentials do not depend solely upon the dilaton field (j> but 
rather upon the so-called shifted dilaton, i.e. <p = <j> — log \/—G s [41, 42, 43] (see also [49]) where 
G s is the determinant of the four-dimensional metric in the string frame and in the cosmic time 
coordinate. 

Therefore the starting point of the present analysis is a generally covariant action supplemented 
by a non-local dilaton potential [46, 47, 48]. The problem will be directly analyzed in the Einstein 
frame metric where the dilaton field and the Einstein-Hilbert term are decoupled and the relevant 
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action is [47] 



where 



S 



d 4 x^—g 



R 



+ 



Ul ' 2 



\g ap d a ct>dpct> - W(4>) 



(2.1) 



e ~* = e ^/ 2 f d 4 y(v^9^d^d^) S(<f>(x) - <j>(y)). 



(2.2) 



In the following the prime denotes a derivation of the corresponding function with respect to his 
own argument (so, for instance, V' will denote a derivation of V with respect to and so on). 



The evolution equations derived from the action (2.1) are then [47] 

^ = 4[W) + W,«7)1, 



2 d<t> 



=Xi - v% 



= 0, 



(2.3) 
(2.4) 



where Q£ is the Einstein tensor already introduced in Eq. (1.2). In Eq. (2.3) T^((f>) and f^(<p,g) 
are, respectively, 



T^) = d^d^-U^d c 

f»(<f>,g) = e*(v6l + l^g^d^d^ e^); 
where the induced metric 7 M „ 

has been defined. In Eqs. (2.4) and (2.6) the two integrals 1\ and I2 are: 



(2.5) 
(2.6) 

(2.7) 



Il(x) = ^ / dAy {^ V ') 6 ^ ~ 

Ux) = ^ / ^{^^g^d^d^ s'i^x) - <f>(y)). 



(2.8) 



Equations (2.3) and (2.4) can now be written for the case of a four-dimensional background 
geometry characterized by a conformally flat Friedmann-Robretson- Walker metric and by a homo- 
geneous dilaton field: 

g^ = a 2 {T)r]^, = 0(t). (2.9) 

Using Eq. (2.9), the set of homogeneous equations stemming from Eqs. (2.3) and (2.4) can be 
linearly combined and the relevant set of evolution equations, in units 2^p = 1, becomes 



e <f, a 2 



H =V2 + 6 



■V, 



(2.10) 
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H>- ^ 4- 



a 2 e" 
6 



b" + 2H(j)' = -eV 



v _3_dV 

2 

v _i_dv 
2 



(2.11) 
(2.12) 



where the prime will denote, from now on, a derivation with respect to the conformal time coordi- 
nate r and 7i = a' /a. 

By summing H' (from Eq. (2.11) ) to 27i 2 (from Eq. (2.10)), the right hand side of the obtained 
equation reconstructs, up to a numerical coefficient, the right hand side of Eq. (2.12) so that the 
dependence on the potential terms may be eliminated. The remaining equation becomes then 
particularly simple: 

^K*+0H- < 2 - i3 > 

Direct integration of Eq. (2.13) gives then: 



n + 



(2.14) 



where e is the integration constant. Unfortunately, in the conformal time coordinate the integration 
cannot proceed further with analytical methods. However, Eq. (2.14) can indeed be integrated 
analytically by defining the new time coordinate a, i.e. 

(It 



da 



(2.15) 



Since the relation of r to cosmic time is dt = a(r)dr, the relation of the cosmic time t to the 
a-parametrization is simply dt = a 3 (a)da. Equation (2.15) can then be used into Eq. (2.14) and 
the result is 

where T = (lna) )t7 and its relation to H is simply H = T '/a 2 . By integrating once Eq. (2.16) 
and by naming as a\ the new integration constant the following explicit relation between the scale 
factor and the dilaton field can be obtained 



(2.17) 



Equation (2.17) simply fixes a relationship between a(a) and 4>{a); such a relationship must hold 
for any potential. Conversely, the particular form of the dilaton potential plays a crucial role in 
determining the relative evolution of the scale factor and of the dilaton whose combination appears 
directly in Eq. (2.17). Finally, once the potential is chosen, the Hamiltonian constraint of Eq. 
(2.10) selects the specific value of e. 

To proceed even further it is necessary to write Eqs. (2.10), (2.11) and (2.12) in the parametriza- 
tion defined by Eq. (2.15); in doing so it is useful to notice that 



T 

n z 



n' 



a 2 \ da 



2JP 



(2.18) 
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Consequently, Eqs. (2.10), (2.11) and (2.12) can be written as 



T 2 = 


i 2 e<V 

— H 

12 6 


-V. 


dT 
da 


2 


V- 


IdV 


d^ = 


<t> 6 

— e v a 


V 


IdV 


da 


~ 2a|. 



(2.19) 
(2.20) 
(2.21) 

where £ = d<fi/da. Equations (2.19) and (2.21) stem directly from Eq. (2.10) and (2.12) by using 
Eqs. (2.18). Equation (2.20) follows from Eq. (2.11) by eliminating, in the obtained expression, 
the dependence on T 2 through Eq. (2.19). 

The solutions of the system formed by Eqs. (2.19), (2.20) and (2.21) can now be derived for 
the general class of potentials V = — Vbe a< ^. In this case we have 

e -{a-2)4>/2 



a 6 e<V 



a 6 e 9 



V 



-V a 6 e*+^ 
IdV 



-Vr 



o- 



a 



3(a-2) 



a _ 2 e -(«-2)*/2 

— n — 



a 3(a-2) ' 



(2.22) 



where we used that <f> = <j) — 31na s , a s (r) being the scale factor in the string frame metric. Fur- 
thermore bearing in mind that the scale factor in the Einstein frame is given by a = e^ 2 a s , Eqs. 
(2.22) follow immediately from the definition of the potential in terms of 4>. 

Equations (2.19) and (2.20) can now be combined and the result is 

dT 3 , o a — 2 



On the other hand, Eq. (2.14) allows to express £ in terms of e and JF, i.e. 

- = e - J-. 



(2.23) 



(2.24) 



Eliminating £ from Eq. (2.23) through Eq. (2.24) the whole problem can be reduced to the following 
equation 

AT m - 9 _ 

(2.25) 



da 2 

that can be easily integrated with the result that T is 



T = -- + ^-tanh 



[e 2 - IT 2 - 2eT], 
f3(a + a 2 ) 



(2.26) 



where (3 = y/S e (a — 2) and a 2 is an integration constant. Inserting Eq. (2.26) into Eq. (2.24), the 
explicit fom of £ turns out tobe 

'f3(a + a 2 y 



£ = 3e — y/3 e tanh 



(2.27) 
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Equations (2.26) and (2.27) can be integrated again so that the explicit functional dependence 
of the dilaton and of the scale factor upon the time coordinate a becomes: 

j. /o e(cr — 2(7-1 ) 

a(a) = e'^^e ^ 

4>(o) = (j, h + Sea - — ^ In cosh ^ ± ^ . (2.29) 

The consistency of Eqs. (2.28) and (2.29) with the Hamiltonian constraint (2.19) demands that 
the value of e is fixed as: 

3e 2 = V r e (a-2) * b e - 3 ( Q - 2 ) CT i. (2.30) 

Equation (2.28) can now be inserted into Eq. (2.28) to obtain the asymptotic relation (i.e. for 
a — > ±oo) between the a and r. This relation is logarithmic in the sense that, up to numerical 
factors a ~ — In (— r) (for cr and r going to — oo) and a — lnr (for a and r going to +oo). This 
is particularly appropriate for numerical purposes since, in practice, the cr-parametrization acts as 
a logarithmic scale in the asymptotic regions. On the contrary, the central region of the dynamics 
(i.e. a ~ (72 ) is described by a linear time-scale. 

To describe the features of the obtained solution let us fix, for simplicity, <ji = 02 = 0. In short 
the basic feature of the solution are the following: 

• for a > 2 the solution (2.28) and (2.29) describes a smooth interpolation between a contracting 
regime (i.e. T < 0) for a < and an expanding solution, i.e. T > valid for a > 0; 

• for a < 2 the dynamical sequence is reversed since the scale factor expands for cr < and 
contracts for a > 0, as is can be argued from the specific form of T (see Eq. (2.26)); 

• for a = 2 the solution always contract and the bouncing behaviour is not realized. 

If <7i or <72 (or both) are different from zero the solutions are not centered around a = and the 
same discussion applies for an appropriately translated coordinate system. 

Concerning the third item of the above list, it is interesting to notice that, in the case a = 2, 
the explicit solution reads 

a(a) = e~^ b/2 e"^ 2 , <j>(a) = 3 e a + b , 3e 2 = V , (2.31) 

implying that T = — e/2. In Fig. 1 (left panel) the evolution of £(cr) and T(a) is illustrated as it 

emerges from Eqs. (2.26) and (2.27). The time scale is given in terms of Cb = 1/e, i.e. the typical 

bouncing time. In the left plot, the logarithm 6 (to base 10) of the scale factor is illustrated for 

6 In the present paper the natural logarithm will always be denoted as In while the logarithm to base 10 will be 
denoted as log. 



cosh 



(2.28) 
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Figure 1: In the left plot the evolution of T '/a^ and £/°"b are illustrated for different values of 
the parameter a. In the right plot the logarithm (to base 10) of the scale factor is illustrated for 
different values of a. 

the same range of a. The initial conditions for the solutions can be charcterized in terms of <p or 
more precisely e^/ 2 that measures the initial strength of the gauge coupling. We shall be primarily 
interested in the case when the solution contracts for a < and expands for a > (i.e. a > 2). In 
this situation the pre-big bang initial conditions will be reproduced since e^/ 2 <C 1 for a —* — do. 
Notice that 0b measures the strength of the gauge coupling around the bounce, i.e. for a ~ 0. Of 
particular interest will be the case when </>b — 0, i.e. when strong coupling is reached just at the 
bounce. It is clearly always possible to tune 0b in such a way that the gauge coupling is much 
smaller than 1 but we shall try to avoid this tuning and examine the more realistic situation where 

4>b ^ o. 

The (7-parametrization is particularly useful for the numerical evolution of the fluctuations and 
for the self-consistent treatment of the dynamical back-reaction effects. It is sometimes practical, 
in a complementary perspective, to solve for the mode functions directly in the conformal time 
parametrization. As already remarked, in the conformal time parametrization we cannot rely on 
exact solutions. Therefore the idea will be to map the asymptotic solutions in the a parametrization 
into asymptotic solutions in the r parametrization and then solve numerically Eqs. (2.10), (2.11) 
and (2.12). It is therefore necessary to develop the numerical intuition for the evolution of the 
background in the conformal time parametrization. Of particular technical relevance is the precise 
mapping of the initial conditions from the cr-parametrization to the conformal time coordinate. 

To fix the ideas, consider the case a = 4 with the bounce occurs around o\ = 02 = 0. In this 
case the mapping of the initial conditions between the a and the r parametrization can be derived 




Figure 2: The evolution of the scale factor (left panel) is illustrated for different initial values of 
the gauge coupling. In the right plot the analytical approximations for the evolution of the scale 
factor (dashed lines) are compared with the numerical solution. 

from the following chain of relations 

p -2</>b 

3e 2 = Voe 2 ^ = — (2.32) 

where the first equality is nothing but Eq. (2.30) in the case a = 4 (and a\ = 0). The second 
equality is fixed by solving directly the Hamiltonian constraint of Eq. (2.10) in the asypmptotic 
regions (i.e. r — > ±oo). Equation (2.32) fixes the relation between the bouncing time in the a- 
parametrization, i.e. = l/ e i m terms of the bouncing time in the r-parametrization, i.e. t\. In 
Fig. 2 the numerical integration of Eqs. (2.10), (2.11) and (2.12) is illustrated in terms of the scale 
factor. In particular, for r < —t\ the solution is well approximated by 




e -^V5±I>, (2 . 33 ) 



while for r > —t\ the solution is well approximated by 



~ - _ „-<^/2 /2(v / 3-l) 



a(T)=a +] j— , a + = e-^4/ vv ^_ > , (2.34) 

where 0b is the same constant introduced in Eq. (2.29). In the pre-bounce region the evolution of 
the dilaton can be obtained from Eq. (2.29) by using the asymptotic relation between r and a in 
the limit r — > — oo and a — > — oo. The result is: 

4>( T ) = 4> b - (V3- l)ln2-\/31n^ V ^ - ) -Vsin^-^J. (2.35) 
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In the intermediate region, the minimum of the scale factor corresponds to 

a min = 0.61 e^ b/2 . (2.36) 

The "width" of the bounce is illustrated in Fig. 2 (right plot) and can be estimated analytically as 

- = _2^ e - fc . I± = ^ e -i,. (2.37) 
ti at Ti a\_ 

In Fig. 2 the values of t± correspond to the two vertical arrows around the origin. 

By changing the various parameters it appears that both a 2 H 2 and 4>' 2 /2 are not sensitive to 
the value of </>b but rather they depend on t\. In particular, it turns out that 

2 0.61 (j)' 2 7.5 
H 2 ^-^, (2-38) 

In the case a = 2 the solution can be expressed exactly also in the conformal time parametriza- 
tion. Recalling Eq. (2.15) we obtain that 

a = ~\n(-e)T, a(r) = e^^y 73 ^, 4>{t) = 4> h - 3 In (-er). (2.39) 

Thus, unlike the case a > 3 the solution with a = 2 is singular in the conformal time parametriza- 
tion. 

The a parametrization has one more important feature that has been already implicitly treated 
in Eq. (2.18). While the sign of Ti and T is preserved by the mapping r — > a, the sign of Ti' is, in 
general, different from the sign of T a . The simplest way of realizing this aspect is to notice that, 
from Eqs. (2.15) and (2.18) we can write 

T = a 2 U, ^ = a 4 {H' + 2H 2 ). (2.40) 

Consider now the simple case of a radiation-dominated Universe. In this case a(r) = (r + 2t\)/t\ 
for r > —T\. Consequently, from Eq. (2.40), Ti and T will have the same sign but J- c will be 
always positive while Ti' = — 1/(t + 2ti) 2 will always be negative. The cure for this type of problem 
is very simple. The invariant quantities (like for instance) the Ricci scalar are always the same in 
both coordinate systems. So, for instance, in a radiation dominated Universe the background Ricci 
scalar, i.e. R is always zero both in the r and in the a parametrizations. In fact, 

s =-> 2+ *'>=-^-4 < 2 -«> 

where the second equality follows by using Eq. (2.18). Now, it is easy to show that, for a radiation- 
dominated Universe, the evolution equation obeyed by T is exactly T t a = F 2 . 
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It is appropriate to remark at this point that the bouncing solutions obtained here in the 
framework of the low-energy string effective action differ from the ones derived in the presence of 
spatial curvature and in a general relativistic context. For instance in [50, 51] (see also [52, 53] and 
references therein) spatially closed models driven by a massive scalar were studied. As anticipated, 
in this framework the dynamical sequence of events differs from the one discussed here. In [50, 51] 
the Universe first deflates, then bounces and finally inflates 7 . In spite of this statement, the 
treatment the dynamical back-reaction of relic gravitons in closed bounces can be developed by 
following the same methods discussed here. 



3 Energy-momentum pseudo-tensor (s) 



3.1 General considerations 



Recalling Eq. (1.1), the quadratic action for hij can be obtained by perturbing Eq. (2.1) to second- 
order in the amplitude of the tensor fluctuations of the metric. Some relevant expressions can be 
found in the appendix A. Up to total derivatives, the result can be written as 



6&S = S gw = J d 3 x dr ^d^dph}^ , 



(3.1) 



where g a g = a?(T)r] a p and r] a p is the Minkowski metric. 

As discussed in connection with Eq. (1.1), the two polarizations of the graviton can be chosen 
to be = — h<2 = h® and /i 1 2 = = h®. Furthermore recalling the notations of Eq. (1.5), Eq. 
(3.1) becomes, for a single tensor polarization 



= }-Jd?xdT^§ d a hdphg a(S , 



Using the parametrization defined in Eq. (2.15), Eq. (3.2) can be also written as 



d 3 x da 



/f)h\ 2 



(3.2) 



(3.3) 



The associated canonical momentum, i.e. II = d a h allows to obtain the Hamiltonian of the tensor 
modes, i.e. 



H, 



n 2 + a 4 d m hd m h 



(3.4) 



By promoting the classical fields to quantum mechanical operators we have that 8 



h(x, a) 



2(2vr) 3 / 2 



d 6 k[h^a)e~ lK - x + hUa)e llt -% 



7 This class of closed bounces has a long history. For instance in [54, 55] useful analytical solutions were derived 
and exploited. In connection with this problem see also [56, 57, 58]. 



8 Notice that, within our set of conventions, - = and ft* ^ = fig. 
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n(^) = / d 3 k[tl^a)e-^ + ul(a)e^ ' 



(3.5) 



2(2vr) 3 /2 J ~ ■-■■-fcv-/- ■ - fc , 
In the Heisenberg representation the evolution of the field operators can be written as 

* :(*i)F k (aT, (3.6) 



^» = ^i)F fc (a)+a T _-(ai)F fe (a)*, 



I%(a) = a A :(a i )G^(a)+a^( f T i )G fc (a)*, 

where a^(cJi) and a_g(<7i) annihilate the vacuum for a; — »■ — oo. The mode functions F k and G k 
now obey 



— + O F fe = 0, G fc = — — , 
da L da 



(3.8) 



where $7 = ka 2 (a). 



3.2 Mixing coefficients 

In the short wavelength limit, i.e. k\r\ > 1 the mode function F k (a) has an oscillating behaviour. 
This regime characterizes the asymptotics of the evolution of the evolution of the mode functions. 
Initially, for a — > — oo the approximate solution of Eq. (3.8) is simply given by 



F k (a) 



1 



G k {a) = — l=[F + iG\e- i f ndff . 



(3.9) 



/2SI V2Q 

In the limit a — > +oo the solution can be expressed by means of two mixing coefficients c+(fc) and 
c_(fc), i.e. 

1 



G k (a) 



'2SI I 
1 



c + (fc) e - i / nd<T + c_(;fe)e < / nd<r 

(^ + in)c + (k)e- i I nda + (ifi - F)c-{ky$ 



Qda 



(3.10) 



In the long wavelength limit, i.e. k\r\ < 1 the approximate solution of Eq. (3.8) across the bounce 
is even simpler and it is given by 

F k = A k + B k a, G k = B k . (3.11) 
It is relevant to remark here that the condition k\r\ ~ 1 reads, in the u-parametrization, 

nda\ ~ 1. (3.12) 



Due to the presence of the absolute values, there will be two distinct moments, in the dynamical 
evolution, where the condition k\r\ = 1 is verified. Before the bounce, i.e. for a < 0, for a given 
mode k the condition —kr(a) ~ 1 will define the time of "exit" in the a parametrization. The 
notation r(cr) simply means that the specific relation between r and a has to be derived from Eq. 
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(2.15) by using the appropriate expression of the scale factor. After the bounce, i.e. for a > 0, the 
condition kr(a) ~ 1 will define the time of "reentry" in the cr-parametrization. 

Defining conventionally the typical time scale of the bounce as Ob = 1 /e, consider the expression 
of the scale factor given in Eq. (2.28). In units of bouncing times the exit and reentry occur, 
approximately, for 

Cex 1 



Ob y/3 + 1 

0"re 1 , 
= — = In 

Ob y/3-1 



l n [^±l e 0b 2 2/(a-2) 
y/3- 1 



,4>b 9 2 /(«- 2 ) 



(3.13) 
(3.14) 



K 

where we defined k = ka\>, i.e. the wave-number in units of inverse bounce time. Equations (3.13) 
and (3.14) can be easily derived from Eq. (3.12) by using the explicit expressions of a(a) given in 
Eq. (2.28) evaluated, respectively, in the limits a < —a^ and a > ov Equations (3.13) and (3.14) 
are derived by fixing o\ = o<i = in Eq. (2.28) and by using Eq. (2.15). Moreover, because of the 
considerations elaborated in the previous Section, the expressions given in Eqs. (3.13) and (3.14) 
hold for a > 2. 

For modes k < 1, the three solutions given in Eqs. (3.9), (3.10) and (3.11) can be appropriately 
matched in such a way that the field operators are continuous function of a. This procedure leads 
to the determination of the coefficients A k , B k 

A k = ^±=e- ikT °*[l + (F ex + in ex )a ex ], B k = - —L^e^™ (F ex + iU ex ), (3.15) 
V zsi ex V 2iZ ex 

and, ultimately, of the mixing coefficients c±(k): 

= ~ 9 7^= g" ifc(Tre+rex) [(^e " ^ex) 
L y a ^cx^ ^rc 

+i(Q re - tt ex ) - (F rc + ^rcX-Fcx + itt ex )(a Te - a ex )], (3.16) 

C + (k) = - _L^ e -^(^-r re)[(i y ex _ j:^ 
L y a ^cx^ ^re 

+i(O rc + O cx ) + (T xc - iVL xc )(T cyL + i^ cx )(cr rc - cr ex )]. (3-17) 

In Eqs. (3.15), (3.16) and (3.17), the obvious notations f2 ex = tt(a ex ), J- Cx = J-{a cx ) (and similarly 
for Q re and J- re ) have been employed. Equations (3.16) and (3.17) imply that |c+(fc)| 2 — |c_(fc)| 2 = 1 
which is a simple consequence of the unitary evolution. For k > 1 the mixing coefficient is exponen- 
tially suppressed. To assess accurately the form of exponential suppression, a numerical treatment 
is necessary. The strategy here is to integrate numerically Eqs. (3.8) with initial conditions given 
by Eq. (3.9) and to obtain, numerically the mixing coefficients, i.e. 

|c + (A0| 2 -MA0| 2 = i[F^(a)G k (a) - F k (a)G* k (a)}, (3.18) 
|c + (A;)| 2 + Mfc)| 2 = ±{\G k (a)\ 2 + (T 2 + Q 2 )\F k (a)\ 2 

+ F[F£(v)G k (v)+F k (v)G* k (v)]}. (3.19) 
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For the numerical consistency of the whole approach, it must be true for each specific numerical 
integration, that the combination reported in Eq. (3.18), corresponding to the Wronskian of the 
solution, is equal to 1 (since it is chosen to be equal to one at the initial integration time, see 
Eq. (3.9)). The numerical accuracy of the Wronskian normalization condition is a check of the 
consistency of the whole approach. Along the same lines, the combination in Eq. (3.19) must go 
to a constant for a — > +oo. 

A set of first-order equations valid for the mixing coefficients will now be discussed. Defining 

e± = e ±i J ndcT , (3.20) 

it is possible to find that the mixing coefficients obey 

dc+ 1 d In 9 dc_ 1 d In 17 9 . 

—± = — e 2 ,c_, — = e 2 c+. 3.21) 

da 2 da + da 2 da ~ + y J 

Equation (3.21) was derived, in a related context, in Ref. [16]. Since c+(fc) and C-{k) are two 

complex quantities subjected to the condition |c+(/c)| 2 — |c_(/c)| 2 = 1, they correspond, overall to 

three real variables that can be chosen to be 

n k = \c-(k)\ 2 , 

Pk = c+(k)c*_(k)e 2 _ + c.(k)c* + (k)e 2 + , 

q k = i[c + {k)c*_(k)e 2 _ - c* + (k)c-{k)e 2 + ]. (3.22) 
The quantities n k , p k and q k obey the following set of equations 

^ = + 2n k ) - 2Qq k , (3.24) 
da 

^ = 2tt Pk . (3.25) 
da 

The initial conditions should be such that, initially, c + (k) = 1 and c_ = 0, i.e. 

n k (ai) = 0, Pfe(cTi) = 0, q k (<Ti) = 0. (3.26) 

Now the analytical estimates will be corroborated by numerical integration. One of the pleasant 
features of the a- parametrization is that the spectrum of the canonical momenta is strictly constant 
in the long wavelength limit. The spectrum of canonical momenta is the Fourier transform of the 
two-point correlation function of Ii(x, a), i.e. 

, sin kr 



where 



(U(x, a)fl(y,a)) = J dink \5 u (k,a)f^, (3.27) 
\Su(k,a)\ 2 = ^\G k (a)\ 2 . (3.28) 
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Figure 3: In the plot at the left, the numerical result of the spectrum of the canonical momenta 
(thick lines), 5u is compared with the analytical expectation (thin lines) for different values of the 
wave-numbers. In the right plot, the numerical evaluation of the mixing coefficients reported for 
k <C 1 ( see Eqs. (3.18) and (3.19)). In both plots the logarithm (to base 10) of the corresponding 
function is illustrated as a function of a in units of bouncing time crt>- The background used for 
this integration corresponds to the case a = 4, ^ = and a\ = C2 = 0. 

Equations (3.27) and (3.28) can be obtained by using Eqs. (3.5) together with Eqs. (3.6) and 
(3.7). When performing the appropriate quantum mechanical expectation values (in the Heisenberg 
representation) it should be recalled that {a^a}^) = 5^ 3 \k — p). 

Consequently, from Eq. (3.11) and from the second relation reported in Eq. (3.15) we get 



where the quantities at the right hand side are fully determined analytically by using Eqs. (3.13) 
and (3.15) and by recalling Eq. (2.26). In Fig. 3 (left panel), the analytical expectation given 
by Eq. (3.29) is compared with the numerical result obtained by integrating Eqs. (3.8) with 
initial conditions dictated by Eq. (3.10). In Fig. 3 (plot at the right) the numerical result for 
the integration of Eqs. (3.8) with initial conditions dictated by Eq. (3.10) is reported in terms 
of the mixing coefficients. We checked that, within the accuracy of the algorithm, the system of 
Eqs. (3.23), (3.24) and (3.25) with initial conditions given by Eq. (3.26) reproduces the same 
asymptotic results for the mixing coefficients. From the analytical estimate the spectrum of the 
mixing coefficients, in the limit a 3> a^, leads to 



that agrees with the explicit numerical calculation of the same quantity reported in Fig. 3. This 




(3.29) 



c + (A0| 2 + Mfc)| 2 ^ + ^ln 2 . 



(3.30) 
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Figure 4: The suppression of the mixing coefficients is illustrated for the regime k > 1. It is clear 
from both plots that when k increases from k < 1 to k > 1, |c_(/c)| 2 — > 0. The background 
parameters are the same as in Fig. 3. 

statement can be explicitly verified by inserting different values of k and by comparing the obtained 
results with numerical values illustrated in the right plot of Fig. 3. For instance, for k = 10 -2 we get, 
from Eq. (3.30), log [\c + {k)\ 2 + |c_(fc)| 2 ] ~ 1.55, for k ~ 1(T 5 we get log [\c + (k)\ 2 + |c_(/c)| 2 ] ~ 2.13 
and so on. 

While the right plot of Fig. 3 nicely illustrates the case k< 1, Fig. 4 deals with the case k ~ 1. 
In this regime, the mixing coefficient c-(k) is known to be exponentially suppressed. It is clear 
that from the plot at the right that as k approaches 1 the mixing is suppressed in such a way that 
for k ~ 1 we have that (|c + (£;)| 2 + |c_(£;)| 2 ) ~ (|c + (A;)| 2 — |c_(A:)| 2 ) implying that |c_(fc)| 2 ~ 0. 

A relevant implication of the results reported so far is that since the spectrum of the mixing 
coefficients is logarithmically increasing with the wave-number (for k = ka^ < 1), the spectra of 
the field operators i.e. 5^ — /c 3//2 |i ? fc(o")| (related to the Fourier transform of the two-point function 
at equal times) will be violet (up to logarithmic corrections). By parametrizing the spectrum as 
$h — /c^* -1 )/ 2 we have, in the example discussed so far, m = 4 (up to logarithmic corrections). 
Moreover, since the spectrum of the mixing coefficients increases with frequency it is plausible 
to expect that the energy spectrum will be even steeper as a function of the wave-number k. 
Consequently, the energy spectrum will be convergent in the infrared. In the ultraviolet, the 
maximally amplified k, i.e. k ~ 1 will give a natural ultraviolet cut-off of the energy density when 
the (initial) zero-point energy is appropriately subtracted. 
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3.3 Energy and pressure densities of relic gravitons 



Different possibilities for assigning the energy and pressure densities of the relic gravitons will now 
be examined. The first strategy is to compute the second-order corrections to the Einstein tensor 
G£ that appears in Eq. (2.3). The details of the calculation are reported in the appendix and the 
main result for the components of the energy-momentum pseudo-tensor can be expressed, in the 
conformal time parametrization, as follows: 



nti u h H + hd m h M d m h M + h' M h M ) 



T 



where 



T-- 



7 



-d m h ke d m h 



i n 



a 2 4 \ 6 



h' ke h ke - -d m h kt d m h 



,mik£ 



4 + \d m h H d m - - A dih M d ] h M 



(3.31) 
(3.32) 

(3.33) 
(3.34) 



with T,\ = 0. From Eqs. (3.31) and (3.33) the components of the energy and pressure density can 



be easily obtained since, by definition, p^ = Tq and 



Pgw 



-r/3. 



The components of the energy-momentum pseudo-tensor given in Eqs. (3.31) and (3.32) are 
not covariantly conserved. However, since the Bianchi identity V^Gjj = should be valid to all 
orders, we will also have that: 

(3.35) 



^ (2) (v M ^) = o, 



whose explicit form is 



(3.36) 



Recalling now the components of the energy-momentum pseudo-tensor and the results for the 
fluctuations of the Christoffel symbols (i.e. Eqs. (A. 3) of the appendix) we have 



dr 



+ 3H(p gw + p gw ) - 



2(H 2 -H>) (2) k _ 
aHl h fe0_U ' 



that can also be written as 



where 



dp. 



dr 



w + 3H(p gw +V gw ) = 



gw 



(h 2 - n') 
ma? 



h' ke h ki . 



(3.37) 

(3.38) 
(3.39) 
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Equations (3.31), (3.33) and (3.39) imply that the operators corresponding to the energy and 
pressure densities are, in the u-parametrization, 

' \±F[hU. + flh] + a A d m hd m h + fl 2 



Pgw 
Pgw 

Pgw = Pgw + 



1 

3^6 



-5IT 2 + 7a 4 (d m h){d m h) 



(3.40) 
(3.41) 
(3.42) 



where h and II are the canonical field operators defined in Eq. (3.5). The averaged components 
of the energy-momentum pseudo-tensor can be obtained by taking the expectation values of the 
operators defined in Eqs. (3.40), (3.41) and (3.42) and by recalling that the initial state is the one 
annihilated by the creation and destruction operators of Eqs. (3.6) and (3.7). The result is then 
expressed in terms of the appropriate mode functions, i.e. 

d 3 k 



(/}gw) = h I Wf {^[F^GtW + F *MG k {a)\ + Q 2 \F k (a)\ 2 + \G k (a) 

1 f d 3 k 
3^ 



(Pgw) 
CPgw) 



^ / T^^m^l 2 - 5\G k (a) 



(2vr 



(Pgw) + ^ 



d 3 k 
(27p 



dhuF 
da 



[F k (a)Gl(a) + FZ(a)G k (a)}. 



(3.43) 
(3.44) 
(3.45) 



It is relevant to remark here that Eq. (3.8) implies that d a (p gw ) + 3^ r ( / o gw + "Pgw) = 0. The initial 
conditions for the field operators (3.6) and (3.7) for a x — > — oo imply that the contribution of the 
zero-point energy is given by J(k 4 )/(2ir 2 a 4 )dlnk (where we recalled that 0, = ka 2 ). To discard 
this quantity by appropriate subtraction (both in the energy and pressure densities) amounts to 
neglect the contribution of the zero-point oscillations of the vacuum. As already mentioned, there 
are different ansatz for the energy-momentum pseudo-tensor that have been proposed in order to 
treat back-reaction effects of the relic gravitons. It is therefore appropriate to comment here about 
these possiblities. By looking at the form of Eq. (3.2), the authors of Ref. [39] (see also [40]) 
argued that a natural ansatz for the energy and pressure densities of the relic gravitons is the one 
we can derive from the energy-momentum tensor of a minimally coupled scalar field for each of the 
two tensor polarizations. This consideration implies that the energy and pressure densities can be 
written as 



1 



Pgw = ^[U 2 + a A d m hd m h], 



a 
1 



IT 



—d m hd m h 

o 



The averaged energy and pressure densities become then: 



(Pgw) — g 



d 3 k 
(27p 



[\G k (a)\ 2 + n 2 \F k (a) 



(3.46) 
(3.47) 

(3.48) 
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d 3 k 



1 r d 6 i 
cfi J (27 



\Gk(a)\ 2 - ^-\F k (a)\ 2 



Equations (3.48) and (3.49) imply d a (pg W } + 3J^(/5 g 
of -Ffc and G k are the ones given in Eq. (3.8). 



(3.49) 

provided the evolution equations 



In order to keep track of the possible ambiguities related with the definitions of the energy- 
momentum pseudo-tensor we will first perform the calculations within the approach defined by 
Eqs. (3.43), (3.44) and (3.45). Then, for each dynamical quantity, the results will be compared 
with the answers obtained using the approach defined by Eqs. (3.48) and (3.49). 



4 Effective barotropic indices 

As a consequence of the dynamical evolution of the background the energy and pressure densities 
have a steep (violet) spectrum. Using the results of the previous section, it appears, for instance 
that for k < fc max we will have, qualitatively, that the logarithmic spectrum of the energy density 



IS a 



where the initial zero-point energy has been already subtracted. For k > /e max the mixing coefficient 
is not enhanced logarithmically, but rather exponentially suppressed (see Fig. 4 where this aspect is 
numerically illustrated). Consequently it is plausible that the integrals appearing in the averaged 
energies and pressure densities are dominated by the modes that are maximally amplified, i.e. 
k ~ K max . Since K max ~ 1, /c max = l/o^ (see Fig. 4). The direct calculation supports this view 
since integrals of the energy density can be performed numerically. Let us then define, for practical 
reasons, 



/^max 
Ai >2 (.k,a)dlnk, (4.2) 



where the indices 1 and 2 refer to the two different parametrizations of the energy density of the 

relic gravitons discussed respectively, in Eqs. (3.43) and (3.48). In Fig. 5 (left plot) the numerical 

evaluation of p gw is reported. The initial conditions for the mode functions are the ones specified 

in the previous section for eack /c-mode. With the dashed line the energy density discussed in 

Eq. (3.48) is illustrated. With the full line the energy density of Eq. (3.43) is instead reported. 

The integral appearing in Eq. (4.2) is indeed dominated by the upper limit. This aspect is also 

illustrated, in a related perspective from the right plot of Fig. 5 where the contribution of the 

integrand for k <C 1 is reported for a sample of /c-modes. In practice, it is useful to specify a given 

ftmin, i.e. the lower limit of integration appearing in Eq. (4.2). For instance, typical values chosen 

9 In the following sections, for sake of simplicity, the averages (/3 gw ) will be denoted simply by p gw and similarly 
for the the pressure densities. 
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Figure 5: Numerical evaluation (left plot) of the integrated energy density for different values of a. 
The contribution of a sample of fc-modes for the case a = 4. Initial conditions for the background 
are fixed in such a way that o\ = 02 = 0, while the gauge coupling at the bounce is of order 1. The 
initial conditions for the tensor modes are the ones dictated by quantum mechanics (see Section 3). 

for the present examples are k = 10 -12 or k = 10~ 10 . Since the spectrum is dominated by K max , 
different choices K m i n produce the same results for the integrated energy density. According to Eq. 
(3.13), K m ; n determines a minimal cr ex , i.e. o" ex (K m in)- This implies that, to be consistent, the initial 
conditions for the field operators have been given for a\ <C o" ex («; m i n ) so that, at the initial time, all 
the tensor modes were of short wavelength. 

The right panel of Fig. 5 has some interesting features that will now be discussed. Within the 
region limited by the two diagonal lines the corresponding modes satisfy \kt\ < 1. Bearing in mind 
Eq. (4.2), we then have that for modes of long wavelength a 6 p gw is constant. Outside the region 
marked by the two diagonal lines the corresponding modes satisfy \kt\ > 1. In this region o 6 p gw 
increases as a 2 , i.e. /? gw ~ a~ 4 . 

These considerations suggest that in the short wavelength limit the effective barotropic index 
is the one of radiation, i.e. 1/3. In fact, in the short wavelength limit the approximate evolution 
of the mode functions is simply given by Eqs. (3.9) and (3.9). Inserting these solutions either in 
Eqs. (3.43)— (3.45) or in Eqs. (3.48)-(3.49) the results are the same for both parametrizations. In 
the long wavelength limit, according to Fig. 5 we should instead have that the effective barotropic 
index is close to 1. 

To corroborate the numerical result, we can solve, directly in the conformal time parametriza- 
tion, the evolution equations of the mode functions in the approximate background of Eqs. (2.33) 
and (2.35). The (averaged) energies and pressure densities will then be computed and the result 
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expanded in the limit of kr <C 1. 

Following this procedure (see appendix B and, in particular Eqs. (B.29) and (B.30)), from Eqs. 
(3.48) and (3.49), we have for r < -n: 



1 



0_ 

2vr 2 
k 4 



dink 



1 



Pgw 



l r k* 
a* J 2^2 



din A; 



. (— irkr) 
1 

( — 7t/ct) 



+ O(fcr) 
+ O(jfer) 



r < — n 



r < —r\. 



(4.3) 



Thus, from Eqs. (4.3) we have that, for kr <C 1, p gw = p gw . Furthermore, recalling that, for 
r < — Ti, a(r) ~ \J— r, we have as expected that /jgw and scale as a -6 . 

If the energy-momentum pseudo-tensor is derived from the quadratic corrections to the Einstein 
tensor, the ratio between the pressure and the energy densities becomes, for kr <C 1, 



' gw 
Pgw 



= 1 - 



3 [1 - 4 7 + 4 In 2 - 4 In (-kr)] ' 



(4.4) 



where 7 is the Euler-Mascheroni constant. Notice that this result is consistent with the one of Eqs. 
(4.3) with the difference that in Eq. (4.4) logarithmic corrections do appear. 

It is instructive to compare the present case with the one of a sudden transition from a de Sitter 
stage of expansion to a radiation-dominated stage of expansion. The typical time-scale t\ will now 
mark the transition from a de Sitter stage of expansion to a radiation-dominated epoch and the 
analytical form of the scale factor in the two regions is: 



a(r) = , 

T 

t + 2ti 



a(r) 



n 



T < -Tl, 

T > — Tl. 



In this case (see appendix B for further details), Eqs. (3.43)-(3.45) imply, for t < — ti, 

k 4 



Pgw 



a 4 I 



P. 



gw 



2vr 2 
1 

" 3a 4 



din A; 
2vr 2 



1 



7 

2^ 

din A; 



1 + 



7 

2^2 



(4.5) 
(4.6) 

(4.7) 
(4.8) 



where x = kr. The energy-density is positive for kr S> 1 but it becomes negative when the relevant 
modes become larger than the Hubble radius, i.e. for A;t <C 1. The equation of state has an effective 
barotropic index 1/3, for kr S> 1. In the opposite limit, i.e. A;t <C 1 the effective barotropic index 
becomes— 1/3. From Eqs. (4.7) and (4.8) the energy and pressure densities scale as a~ 2 and this is 
consistent with Eq. (3.38) for the equation of state V gw = —p gw /3. 

Consider then the case of modes kr <C 1 for t > — ti i.e. during the radiation-dominated 
epoch. The exact results are reported in the appendix and here only the long wavelength limit 
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will be given. As as the approach defined by Eqs. (3.33)-(3.34) is concerned, the relevant results 
are reported in Eqs. (B.13), (B.15) and (B.16) of the appendix. To lowest order in \kr\ <C 1 and 
| fcn | <C 1 Eq. (B.13) leads to 

5 1 r dink 
6^2 J ^2^4 

In the same limit, from Eqs. (B.15) the following result 



[ dink 2 



7 1 f d In k , 9 , 

**° = %*lw%*> (410) 

can be obtained. We can immediately notice that p gw /Pgw = — f , which agrees with the calculation 
of [36, 37]. The fact that the averaged energy density can have negative values for a limited amount 
of time was noticed, in a related context, in [59]. 

Finally, from Eq. (B.16), to lowest order in \x\ <C 1 and \x\\ <C 1, 



8 1 f dink 2 



Therefore, putting together the results expressed in Eqs. (4.9), (4.10) and (4.11), the following 
chain of equalities holds 

7 

Pgw = -^p g w, (4.12) 

'Pgw = ~^Pgw + rj^Pgw = ~~qPs w ' (4-13) 
These results coincide with the ones obtained in [36, 37]. 

Instead of using the energy-momentum pseudo-tensor derived from the quadratic corrections to 
the Einstein tensor, we could use the energy and pressure densities discussed in Eqs. (3.46) and 
(3.47). The relevant results, in this case, are reported in Eqs. (B.18) and (B.19) of the appendix. 
Again, to lowest order in |x| <C 1 and \x\\ <^ 1, Eqs. (B.18) and (B.19) lead to 

1 r 1.4: 2 

' ' " ' " ' - ■' (4.14) 



i r k A „ x 2 



1 6 2 



(4.15) 



where the ellipses stand for higher-order corrections. While in the case of Eq. (4.9) the averaged 
energy density is negative, in the case of Eq. (4.14) the averaged energy density is positive. More- 
over, from Eqs. (4.14) and (4.15) we have that the effective equation of state in the long wavelength 
limit is given by 

Pgw = -gPgw, (4.16) 

which is the same of Eq. (4.13). 
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Therefore, we can conclude that both in the bouncing case and in the transition from de 
Sitter to radiation the two parametrizations of the energy-momentum pseudo-tensor give the same 
information. In particular, Vgw defined by means of the corrections to the Einstein tensor obeys 
the same equation of state obeyed, in the approach of Eqs. (3.46) and (3.47) by p gw . 



5 Iterative calculations of back-reaction effects 



The evolution equations including the dynamical effects of the produced gravitons reduce to an 
integro-differential system whose numerical solution may be obtained by means of iterative methods 
that we are now going to describe and exploit. One of the important qualitative features of the 
bounce solutions discussed in Section 2 was expressed by the Eq. (2.13) whose form becomes now 



d ( r n a\ 



(5.1) 



where with pg„ and p gw we denote the averages of the components of the energy-momentum pseudo- 
tensor of the relic gravitons. To this equation the evolution equation for £ should be added, i.e. 



da 



V — 



dV_ 



(5.2) 



which is (formally) independent on the energy and pressure densities of the created gravitons. The 
dynamical quantities appearing in Eq. (5.1) and (5.2) are subjected to the following (generalized) 
Hamiltonian constraint: 

-V, (5.3) 



T "^ + 12 + 



6 



that has to be satisfied along the different steps of the numerical solution. Recalling the examples 
given in Fig. 5, p gw and p gw are computed by integrating the evolution equations of the mode 
functions, i.e. Eqs. (3.8). In this sense, Eq. (5.1) is actually an integro-differential equation whose 
explicit form is indeed 

B(Q,a) d\nk, (5.4) 



d_ 

da 



where the quantity B(£l,a) depends on the specific form of the energy-momentum pseudo-tensor. 
In particular, subtracting Eq. (3.49) from Eq. (3.48) 



2k 3 

B(n,a) = —[n 2 (a)\F k (a)\ 2 - Q(a)}. 

Similarly, taking the difference between Eq. (3.43) and (3.44) 
2k 3 



(5.5) 



B(n,a) 



3vr 2 



{2\G k (a)\ z - &{a)\F k {a)\ z + ^[F k {aYG k {a) + F k {a)G k {a)* - fi(<r)]}. (5.6) 
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In both equations Fk(a) and Gfc(cr) satisfy Eqs. (3.8). Note that in Eqs. (5.5) and (5.6) the 
zero-point energy and pressure densities have been subtracted as discussed in Section 4. 

Equation (5.4) can be solved by iteration 10 . In fact, from Eq. (5.4) the following hierarchy of 
equations can be deduced: 

*b(*) + ^ = 0, (5-7) 
jri ( <J )+ 2 = l7 da> Jo Mn{a'),<j')d\nk, (5.8) 
T 2 (a) + ^ = \ f da' jf**"" Bi a') d\n k, (5.9) 

and so on. 

The subscripts appearing in the various dynamical quantities of Eqs. (5.7), (5.8) and (5.9) 
denote the order of the iteration. The source term appearing in the hierarchy is computed from the 
previous order. So, to zeroth order, i.e. Eq. (5.7), To and £o are simply the solutions discussed in 
Section 2. These solutions will be used to obtain numerically the mode functions that will determine 
23o(f2, cr). To first order, T\ and £i will be determined from the the zeroth order source, i.e. Bo(£l, a). 
This will lead to the first-order mode functions allowing the quantitative determination of £>i(f2, <r) 
and so on. 

Based on the results of the previous section, the strongest effect is expected from modes of high 
wave- number (i.e. small wavelength). The rationale for this expectation is that in the long wave- 
length limit a 6 (p gw — p gw ), is smaller than the energy density of the dilaton and also approximately 
vanishing. In the short wavelength limit a 6 (p gw —p gw ) increases as a 2 and it is expected to dominate 
the background at sufficiently late times: as time goes by more and more modes will contribute to 
the ultraviolet branch of the energy spectrum. For the numerical solution of the problem the (stiff) 
Rosenbrok method [61] has been used. In Fig. 6 the back-reaction effects are reported for the case 
of a zeroth-order solution with a = 4. With the dashed line the zeroth-order iterative solution is 
illustrated. With the full line the combined effect of the first-order iteration is illustrated. In a 
radiation-dominated Universe the evolution equation obeyed by F(o) is simply given by 

(5,0) 

where p r is a generic radiation fluid with p r = p r /3. As already noticed in connection with Eq. 

(2.41) the first of the two relations appearing in Eq. (5.10) implies that the Ricci scalar vanishes 

in a radiation-dominated Universe. The behaviour of the solution for a 3> <7t>, illustrated by the 

full lines in Fig. 6, corresponds, indeed to the one of a radiation-dominated Universe in the a 

10 It should be mentioned that an iterative approach for the solution of the back-reaction problems has been also 
invoked by the authors of Ref . [60] . 
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Figure 6: The evolution of the scale factor (left plot) and of T (right plot) in the case a = 4. 
The back-reaction effects are implemented using the Eq.(5.5). The parameters of the zeroth order 
solution are fixed by requiring o\ = 02 = and </>b = 0. 

parametrization. In this parametrization, in fact T ~ (a — cr*) -1 and T '/a is constant. This non- 
standard evolution of the scale factor becomes more familiar by translating the result from the 
<T-parametrization to the r parametrization. Recalling, in fact, that a(a) 2 da = dr, we do find that 
r ~ (a — o"*) -1 . This occurrence implies, as expected, that a(r) ~ r which is the well known result 
for a radiation-dominated Universe in the conformal time parametrization. 

In Fig. 6 the back-reaction effects have been implemented by using Eq. (5.5) which comes from 
Eqs. (3.48) and (3.49). In Fig. 7 the numerical calculation has been instead performed by using 
Eq. (5.6). By comparing Fig. 6 with Fig. 7 the main visible effect is an enhanced oscillating 
behaviour of the first-order iteration after the bounce. However, the evolution of the scale factor 
and of J- is clearly the same. 

The iterative method described in the present Section can be also applied in the conformal 
time parametrization. We will omit here the technical details but the general idea is to start with 
the asymptotic solution valid for r -C —t\ (which is, for instance the one provided by Eqs. (2.33) 
and (2.35) ) and then integrate numerically both for the background and for the fluctuations. The 
main technical difference between the integration in the a and in the r parametrization is the 
following. In the a parametrization the analytical solution of Eqs. (2.26) and (2.27) can be used 
in the zeroth-order iteration. In the r parametrization already the zeroth order iteration requires 
numerical treatment. Moreover, in the r parametrization the evolution equations for the mode 
functions are formally different and are the ones reported in Eq. (B.22) of the Appendix B. 

In Fig. 8 the integration of the iterative problem is illustrated in the r parametrization. With 
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Figure 7: The same quantities reported in Fig. 6 but using implementing the back-reaction effects 
through the energy and pressure densities stemming from the quadratic corrections to the Einstein 
tensor (i.e. Eq. (5.6)). 

the dashed line the zeroth-order solution is reported. With the full line the result of the first 
iteration is illustrated. In the left plot of Fig. 8 the parameters are chosen to be the same ones 
of Figs. 6 and 7. In the right plot the value of the gauge coupling at the bounce is chosen to be 
much smaller than 1. This implies that the scale factor after the bounce is much larger than in 
the case when the gauge coupling is of order 1 at the bounce. This can be understood analytically 
since the asymptotic zeroth-order solution (i.e. Eq. (2.34) leads to a scale factor whose absolute 
normalization increases as the gauge coupling at the bounce diminishes. It is evident from Fig. 
8 that the first-order iteration already leads to a scale factor that increases linearly in conformal 
time. For both plots in Fig. 8 the back-reaction effects has been parametrized through the analog 
of Eq. (5.5), appropriately translated in the r parametrization. 

6 Concluding remarks 

Dynamical back-reaction of relic gravitons is relevant in different contexts, so that the nature of 
the induced physical effects may vary. A particularly interesting situation is the one of bouncing 
solutions of pre-big bang type where back-reaction effects are known to be important already from 
qualitative estimates. One of the purposes of the present investigation is to propose a theoreti- 
cal scheme where more quantitative predictions could be derived. The present findings are then 
useful for implementing a smooth exit from the pre-big bang phase with consequent production 
of radiation. The bouncing solutions examined here may also be viewed as a practical theoretical 
laboratory where the methods for the analysis of the back-reaction effects of relic gravitons may 
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Figure 8: Integration of the iterative problem in the conformal time parametrization. In the left plot 
the parameters of the zeroth-order solution coincide (in the limit r — ► — do) with the parameters 
chosen for the integrations reported in Figs. 6 and 7 (for a — ► — do). 

be tested. In this respect, the reported results suggest that the role of the produced gravitons is 
crucial for determining the asymptotic state of the solution at late times. 

The strategy followed in the present analysis relies on the interplay between analytical methods 
and numerical calculations. We have been interested in situations where the Universe undergoes an 
accelerated contraction before the bounce, while, after the bounce, the zeroth-order solution exhibits 
a decelerated expansion. After quantizing the tensor modes of the geometry, the initial values of the 
field operators have been set well before the bounce. At that time all the gravitons were of short- 
wavelengths in comparison with the typical time scale of the background evolution. Then, the field 
operators have been evolved in the Heisenberg representation and the mixing coefficients have been 
accurately determined both analytically and numerically. This step allowed to understand what 
kind of regularization is required for the energy and pressure densities of the produced gravitons. 
The second step has been to adopt a consistent ansatz for the gravitational energy-momentum 
pseudo-tensor. Two complementary definitions of the energy and pressure densities of the relic 
gravitons have then been scrutinized and compared both analytically and numerically. On the 
basis of the reported results, these two different approaches lead to compatible quantitative results. 

The self-consistent analysis of back-reaction effects leads naturally to an integro-differential 
problem. Consequently, an iterative method for its solution has been proposed and tested in the 
specific dynamical framework of bouncing solutions that arise in the context of the low-energy string 
effective action supplemented by a non-local dilaton potential. It turns out that the treatment both 
of the background and of the tensor fluctuations becomes rather simple by selecting a different 
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parametrization of the time coordinate which does not coincide with the usual conformal or cosmic 
times. In this time parametrization analytical bouncing solutions can be derived in the Einstein 
frame metric. 

By looking at the structure of the low-energy string effective action, it is possible to infer that 
not only relic gravitons but also other massless fields may lead to similar effects whose analysis is 
certainly one of the possible developments of future studies. In this respect, the present investi- 
gation dealt with the minimal field content of the model with the aim of achieving a reasonably 
accurate description of the underlying physical processes. Thus, some of the ideas discussed here 
will hopefully be useful also for the description of dynamical systems with a richer field content. 
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A Energy-momentum pseudo-tensor 



From Eq. (1.1) we have 

S^gij = -a 2 hi,, hi = dih) = 0, (A.l) 

where, as in the bulk of the paper, the subscript refers to the tensor nature of the fluctuation while 
the superscript denotes the perturbative order. Equation (A.l) implies that 

4V = %, «SV - (a.2) 

or a z 

Consequently, the fluctuations of the Christoffel connections to first and second order become: 





^(h' ij + 2Hh ij ), 


^rl = 


V 


«t (1) r&- = 


I(^ + 5A fe -5%), 


*I 2 H = 






-h %e [deh jk - d k hj£ - djh ke ] 



(A.3) 

where, as in the bulk of the paper, the prime denotes a derivation with respect to the conformal 
time coordinate. The fluctuations of the various components of the Ricci tensor to first and second 
order are then: 

S^Rij = \W J +2Hh> lj -V 2 h lj } + (H > + 2H 2 )h ij , (A.4) 
<5 t (2) i?oo = ^/-yM«' + ^V\, (A.5) 
S^Rij = ^h M [d k d e hij - d k djh ei - d k dih je ] 



2 

1 1~L 

-dj[h ke (d e h ik - d k h u - dih M )} - —h kl h' kl 8 i3 



+ -^h)h' ti + -h\h' tj - -h k h' ik - -h) h ik - -h k h' kj - -h k h kj 

- - [dih{ + d k h\- d e h lk ] [d e h k + djh k t - 8 k h je ] . (A.6) 



The Ricci scalar is zero to first order in the tensor fluctuations, i.e. 5^ R = 0. This is due to the 
traceless nature of these fluctuations. To second-order, however, S^R / and its form is: 

+ ^{-ld t [h k \d e hi-d k h}-d*h ke )} 

- x - [fyhidehl - dih{d k h u + d k h u d e h k - tfhik&ht + d% k d k h u \ } . (A.7) 
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Recalling the form of the Einstein tensor (see Eq. (1.2) ), 

?( 2 )rr„„ - _ r(2) 



we obtain 



4r 00 = + J(/t' M ^' + dihud'h^) + Poo, 



where 2?oo is a total derivative, i.e. 



V Q0 = ^d i [d i h k %-2d k hlh hi ] 



(A.8) 
(A.9) 

(A.10) 



From Eqs. (A. 6) and (A. 7) it is also possible to write: 



s[ 2 ^ Rij 



-(d k hid k h je + a l h ik a e h k j ) 



7 a i h ke d j h ke 



h k %t5 ij + ^-(h% + h<ih' ej ) 



H \k> ,,! ,s /! 1,1 1,1 1,1 ; 



H (h k 'h lk + h k 'h kl ) - - A (h k 'h' lk + h k 'h>) + Vi 



si 2) R 



l 



^h' ke h ke ' + Hti ke h k 



\d^ U &h kl 



where T>ij and Vr are further total derivative 



T>ij = ^d k [h ke (dehij - djhu - dih je )] - \^dj\h k \dth ik - d k h u - dih M )} 
1 



- -d e [d k hW + h ik &'hZ]. 



1 



V R = d i [h k ^h H ] + -d e [d i h k %-2d k hlh% 
Therefore, up to total derivatives, the following result holds: 



1 



and 
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4 Too = Hh' kl h + -{h' u h M ' + dih u &h M \ 



%3 = l^[dmh ke d m h ke - h' ke h H '} + l -h k 'h' lk + -^hud^ - \d k h\d k h lr 



(A.ll) 



(A.12) 



(A.13) 



(A.14) 

To pass from doubly covariant indices to mixed ones, it is useful to recall that, to second order, 

8^Ql = 8® \g ua Q m \ = b^g™Q m + g ua 5 {2) Q, a + <5 ( V^% a . (A.15) 



By looking at the form of the specific terms arising in the previous equation it is clear that T§ = 
5 00 7oo and that T? = ~gi k T ki . The expressions for and 1? are the ones reported in Eqs. (3.31) 
and (3.32) (see also Eqs. (3.33) and (3.34)). These expressions coincide with the ones obtained, for 
instance, in [36, 37] and are also consistent with the ones of [34, 35]. 
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B Analytical estimates of the effective barotropic indices 



Since the Hamiltonian defined in Eq. (3.4) is time-dependent, it is always possible to perform 
time-dependent canonical transformations, leading to a different form of the Hamiltonian that will 
be classically equivalent to (3.4). This procedure correspond to drop total time derivatives from 
the corresponding classical action [62]. 

From the action (3.3), defining the new field fx = ah and dropping a total time derivative the 
following Hamiltonian can be obtained in the r-parametrization 



7T 2 + (5,/i) 2 - (a. 2 + n')^ 



(B.l) 



where ir = ji 1 . This form of the Hamiltonian is particularly convenient for the studying the time 
evolution of the field operators in the conformal time parametrization. By appropriate Fourier 
transforms of the field operators (see Eq. (3.5)), the analog of Eqs. (3.6) and (3.7) are 



Afc(r) = kin) Mr) + 11^)^), 
*fc( r ) = k( T 09k(r) + bK(T0gUr), 

where the new mode functions now obey 

a" 1 



k - 



fk — 0, g k — f' k . 



(B.2) 



(B.3) 



In the limit t\ — > — oo, the operators bAj\) and b_AT\) annihilate the same initial vacuum state as 
the one annihilated by o^(ci) and a_g(o"i) in the limit o\ — ► — oo. 

In this framework all the calculations discussed in the cr-parametrization can be reformulated. 
The expectation values of the energy and pressure density can be computed. For instance from 
Eqs. (3.31) and (3.32) we have 11 



1 



(Pgw) — 4 



a 
1 

3a 1 



(X)3 
d 3 k 



{{k 2 - m 2 )\f k {r)\ 2 + \g k {r)\ 2 + m[f k {r)gt{r) + / fc *(r) 5fc (r)]}, (B.4) 



(2vr) 3 

+m[f* k (r)g k (r) + f k (r)gUr)}}. 



{\f k ( T )\ 2 (7k 2 -5H 2 )-5\g k (T)f 



CPgw) — (Pgw) + 



4(H 2 - W) 
3Ha 4 



d 3 k 



[f k (r)g* k (T) + f* k (r)g k (T) 



(B.5) 
(B.6) 



In similar terms, the expectation values of the energy and pressure densities can be obtained within 

the approach of Refs. [39, 40]. 

11 In the bulk of the paper, for practical reasons, we dropped the expectation values for the averaged quantities. In 
this appendix, to make clear the computational procedures, the expectation values will be restored. 
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Consider now the case of a transition from a de Sitter stage of expansion to a radiation- 
dominated stage of expansion, i.e. Eqs. (4.5) and (4.6). For r < — n, 



1 n *' 
— -1 + - 

x z x 



e~ ix , (B.7) 



V2k\ x, 

where x = kr. For r > — t\ the canonical fields can be expressed as 

Afc(r) = S^(r) + 6 f _^(r), %(r) = b % G k {r) + S^-G^r), (B.8) 
where the mode functions are now 

F k (r) = J= e -^ +2 ^\ G k {r) = -Jt e -^+^. (B.9) 



Since the field operators must be continuous we have 



/fc(-ri) = c+(k)F k (-n) + c-(k)F£(-n), 

9k(-ri) = c + (k)G k (-n) + c-(k)G* k (- Tl ). (B.10) 



so that c + (k) and c-(k) can be determined as 

„, e 2ixi (2x 2 - l + 2zxi) 1 .„ , 

c + (k) = U t c-{k) = ^ (B.ll) 

where x\ = kr\. It can be verified immediately that |c + (fc)| 2 — \c-(k)\ 2 = 1. 

It is now instructive to consider the behaviour of the effective barotropic indices in the two 
regimes r < — n and r > — t\. Furthermore, in each regime, we shall be interested in the relative 
weight of the short and long wavelength modes. This analysis is useful for the comparison with 
the case of bouncing models. Since during the de Sitter phase, H' = H 2 , Eqs. (B.5) and (B.6) 
imply that (pgw) = {Pgw}- Furthermore, if we will have that 3(p gw ) = (p gw )- During the radiation- 
dominated phase, i.e. r > — t\, the short wavelength modes also obey a radiative equation of state 
and the energy and pressure densities explicitly depend upon the mixing coefficients, i.e. 

</W> = ^/ |^[2M*)| 2 + 1], 
1 f krfik 

(V gw ) = (p gw ) = —,J ^3 [2|c_(fc)| 2 + 1], (B.12) 

where, for simplicity, the vacuum contribution has not been subtracted. Notice that, now (i.e. in 
the radiation dominated epoch and for kr S> 1) , the equality of {Vgw) and (p gw ), is a consequence 
of the explicit form of the mode functions in this regime, i.e. F k G\ + F k G k = 0. 

In the long wavelength limit, i.e. kr <C 1 and during the de Sitter phase we obtain, instead, the 
same expressions reported in Eqs. (4.7) and (4.8). Consider finally the last case, i.e. long wavelength 
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modes during the radiation dominated epoch. After making use of the explicit expressions of the 
mixing coefficients we will have, for the averaged energy density, 



1 f k ct 
(Pgw) - ~ 4 J 



1 r k d 3 k(l 



(2vr) 3 I 2 
— sin2(x + xi) 



7 



(x + 2x\ 



3(2x? - 1) 



+ 



+ cos2(x + xi) 



2(x + 2xi)xf 2(x + 2xi) 2 xf 
3 7(2x 2 - 1) 



.(x + 2xi)xf 4x|(x + 2xi) 2 . 
Consequently, to lowest order in |x| <C 1 and |xi| <C 1 the averaged energy density will be 

5 1 r dink 



5 1 f dink 2 



(B.13) 



(B.14) 



which is the result reported in Eq. (4.9). Along similar lines the averaged pressure density can be 
obtained: 

1 [k 4 dlnk(l 



(Pgw) 



3a 4 / 



2vr 2 1 2 
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4xT 
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(x + 2x 



(x + 2xi) 2 _ 

-y ^ [(2xf ~ 1) cos 2 ( a; + ^i) + 2a: i sin 2 ( :r + x i)\ 



— —t[(2x\ - l)sin2(x + xi) - 2xi cos 2(x + xi)] 1, (B.15) 
iK J 



CPgw) (Pgw) 
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1 \ (2x\ ~ 1) 



2xf 



2xxf 



cos 2(x + xi) 



sin2(x + xi) 



(x + 2xi)xf 

implying, to lowest order in |x| <C 1 and \x\\ <C 1, 

7 1 f dink 



('Pgw) (Pgw) 



6 a 2 J 2tt 2 t 1 4 ' 
which are the results reported in Eqs. (4.10) and (4.11). 



8 1 f dink , 9 



9 a 2 7 2tt 2 t! 
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(B.16) 



(B.17) 



Instead of using the energy-momentum pseudo-tensor derived from the quadratic corrections to 
the Einstein tensor, we could use the energy and pressure densities proposed in [39, 40] (see also 
[6]). In this second case the expectation values of the energy and pressure density become: 

k 4 



{M = ^ I 2^ dlnk {l 
(2x\ - 1 



2 + 



+ cos 2(x + xi) 



(x + 2xi) 2 
1 



1 + 



2x? 



4x 4 (x + 2xi) 2 xf(x + 2xi) 
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+ sin 2(x + x\) 



(2x 



+ 



2x\{x + 2xi) 2xf(x + 2xi) 
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k' 1 
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(B.18) 
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(B.19) 



If we now expand the previous expressions for |xi| <C 1 and for |x| « 1 we get to the results 
reported in Eqs. (4.14) and (4.15). 

In the following we are going to cross check some of the results obtained in the case of the a 
parametrization but working directly in the conformal time coordinate. It should be clear that 
the specific parametrization is not crucial to obtain a given result. However, it can happen that, 
for practical reasons, a given parametrization has to be preferred since it leads either to more 
transparent or to swifter results. The idea will be here to solve the evolution of the mode functions 
exactly but keeping an approximate form of the background in the post-bounce and in the pre- 
bounce region. In the conformal time parametrization the action given in Eq. (1.4) reads 



S m = \j d 3 xdTa 2 [(d T h) 2 -a 2 {d m h)\ 



(B.20) 



whose associate Hamiltonian is 



H ^ = \j d3x + fl2 ( 9 ™ ( B - 21 ) 
with canonical momentum given as II = a 2 h'. Therefore in this case the mode functions obey 



f' + 2Hf' k + k 2 f k = 0, ~g k = a 2 f k 



(B.22) 



Consider, for instance, the case given in Eqs. (2.33), (2.34) and (2.35). The form of the mode 
functions before the bounce (i.e. for r < — n) will then be 



fk(r) 



e in/4 H^\-x), 



9k 



\/2a_ 

(r) = -a'^k^/^H^i-x), 



(B.23) 



where, as in the previous example, x = kr; H^(z) and H^\z) will denote, in the following, the 
Hankel functions of index v and argument z [63, 64]. For r > ti, recalling Eq. (2.34), the evolution 
of the mode functions can be written as: 



Fk(r) 



V2 a+ 



c+{k) e™l A H§\x) + c-{k) e~™/ A H {2) \x) 



G k (r) = -a 2 (r) 



a + 



c + {k) e iw / 4 H{ 1 \x) + c-(k) e~ i7T / 4 H{ 2 \a 



(B.24) 
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It is now instructive to pretend that the analytical form of the mode functions is unavailable in 
the intermediate region, i.e. —t\ < t < r±. Therefore let us try to excise a thin region around the 
bounce and compute the mixing coefficients by assuming 

/fc(-Ti) = F fc (n), g k (-Ti) = G fc (n). (B.25) 

This allows to determine 

c+{k) = -ijX 1 [H[ 2 \x 1 )H^(x 1 ) + H^(x 1 )H[ 1 \x 1 )}, 

c-(k) = —x^ix^H^ix!), (B.26) 

Notice that from the previous formulas it is possible to obtain, after some trivial algebra, that 
|c + (fc)| 2 - |c_(/c)| 2 = 1. Recalling now that [63, 64] 

fl-J 1) (xi)~t-lnxi, H^ixi) — (B.27) 

7T TTX± 

we find that for x\ <C 1 

c-(k) -In xi (B.28) 

7T 

Now this result is (qualitatively) similar to the one obtained, by a totally different approach, in 
Eqs. (3.16), (3.17) and (3.30). By doing a more careful comparison it is possible to show that the 
method discussed in Section 3 of the present paper is the most accurate, also from a conceptual 
point of view. However, the path followed in the present appendix is able to reproduce the gross 
features of the mixing and, therefore, can be used for semi-quantitative estimates. 

It is interesting to notice that, within this approach, the results of the effective barotropic 
indices (numerically obtained in the text) can be checked. For instance, in the approach of [39, 40] 
the energy and pressure densities before the bounce become, in the present example 

(/}gw) = ^/|^ dln 4l ( " x)[Jl( " a;)2 + yi( " a;)2 + Jo( " x)2 + yo( " a;)2] }' 
iM = ^J^2d\nk^(-x)[J 1 (-x) 2 + Y 1 (-x) 2 } 

-y 2 (-x)[M-xf + y (-*) 2 ]}, (B.29) 

where J v (z) and Y u (z) are the Bessel functions of order v and argument z [63, 64]. In the small 
argument limit, from Eqs. (B.29) we do get: 

(Pgw) ^ (Pgw) = ±f ^2 d ^k(-^), (B.30) 

which is one of the results discussed in Section 4. Both (pgw) and (pgw) scale as a" 6 (since a(r) ~ 
y/—r)as it is consistent with the equation of state (p gw ) — (Pgw)- 
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The same analysis can be preformed when the energy-momentum pseudo-tensor is obtained 
from the quadratic corrections to the Einstein tensor. The result is 

k 4 



(Pgw) = 



1 



2vr 2 



d\nk\ir{J (-x)J 1 (-x) +Y Q (-x)Y 1 (-x)] 



j[H-xf + M-xf + y (-x) 2 + yi(-x) 2 ]}, 

1 f k 4 { 
(Pgw) = ^4 / ^ dlnk \ n[Jo(-x)Ji(-x) +Y (-x)Y 1 (-x)] 

jx[U-x) 2 + Yo(-x) 2 + M-x) 2 + Fi(-x) 2 ]}, 



(Pgw) (Pgw) 



2vr 2 



dlnM7r[J (-x)Ji(-x) +Y (-x)Y 1 (-x)}. 



Let us now expand the previous expressions for |x| <^ 1. The result is: 



(Pgw) 
(Pgw) 
(Pgw) — (Pgw) 



a 4 I 



dink 



— -(l-47 + 41n2-41n(-x)) 

TTX 



a 4 I 



k^_ 
2vr 2 
k 4 5 

k 4 2 

^rdlnk — [2 7 -21n2 + 21n(-x)l. 
2ir z nx 



(B.31) 



(B.32) 



where 7 is the Euler-Mascheroni constant [63, 64]. From Eqs. (B.32) it is possible to obtain that in 
this parametrization of the energy-momentum pseudo-tensor, indeed, logarithmic corrections to the 
effective barotropic index arise, as pointed out in Section 4. From the examples reported here it is 
clear that also the conformal time parametrization can be employed for the self-consistent analysis 
of back-reaction effects. Indeed, in Section 5, explicit numerical integrations have been reported 
also in the r parametrization. It is however true that, in the context of bouncing models, the a 
coordinate allows to obtain the wanted results much faster both analytically and numerically. 
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